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Abstract 

We introduce three distinct, yet equivalent, optimization procedures for the Fourier 
Spectral Method which increase its accuracy. This optimization procedure also allows us 
to uniquely define the error for the cases which are not exactly solvable, and this error 
matches closely its counterpart for the cases which are exactly solvable. Moreover, this 
method is very simple to program, fast, extremely accurate {e.g. an error of order 10"^^" 
is usually obtainable as compared to the exact results), very robust and stable. Most 
importantly, one can obtain the energies and the wave functions of as many of the bound 
states as desired with a single run of the algorithm. We first thoroughly test this method 
against an exactly solvable problem and then apply it to two problems which are not 
exactly solvable, all for finding the bound states of the time-independent Schrodinger 
equation. We present a detailed comparison between the results obtained by this method 
and some of the more routine methods. 



1 Introduction 

Eighty years after the birth of quantum mechanics [1], the Schrodinger's famous equation 
still remains a subject for numerous studies, aiming at extending its field of applications and 
at developing more efficient analytic and approximation methods for obtaining its solutions. 
There has always been a remarkable interest in studying exactly solvable Schrodinger equations. 
At this point, we have to state that traditionally the term "exactly solvable" has been used 
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in a well-defined mathematical sense, meaning that the eigenvalues and eigenfunctions of the 
Hamiltonian under consideration may be expressed in an explicit and closed form [2]. In 
this sense, the exact solubility has been found for only a very limited number of potentials, 
most of them being classified already by Infeld and Hull [3] on the basis of the Schrodinger 
factorization method [4] , which in turn appeared to be a rediscovery of the formalism stated 
nearly 120 years ago by Darboux [5]. However, a vast majority of the problems of physical 
interest do not fall in the above category when we formulate a more or less realistic model for 
them. Then we have to resort to approximation techniques which can be analytic or numeric. 
Examples of approximate analytic techniques would be the usual perturbation method, the 
semiclassical or WKB approximation, and the variational method [6, 7]. In the perturbation 
method the problem is divided into two segments. The main segment is supposed to be 
exactly solvable. The second segment is supposed to modify the solution obtained in the main 
segment only very slightly. This modification can be obtained analytically for any desired 
degree of accuracy. On the other hand, the numeric solutions could be either perturbative in 
nature or completely numeric. Indeed, the relevant Schrodinger equation can always be solved 
numerically, and in view of the immensely increased computational power, this numerical 
solution could be extremely accurate. Such solutions could serve the dual purpose of being a 
substitute for the exact solution and being used as a comparison with experimental results, 
when applicable, whose accuracy are also rapidly on the rise. The need for such methods 
have stimulated development of more sophisticated integration approaches, e.g. embedded 
exponentially-fitted Runge-Kutta [8] and dissipative Numerov-type [9] methods, as well as 
interesting techniques, such as a relaxational approach [10] based on the Henyey algorithm 
[11], an adaptive basis set using a hierarchical finite element method [12], and an approach 
based on microgenetic algorithm [13], which is a variation of a global optimization strategy 
proposed by Holland [14]. To be more specific, Amore [15] has used the variational sine 
collocation method to obtain the bound state solutions for a Schrodinger equation and has 
obtained accuracies of order 10~^^. Also Jentschura [16] has used the Instanton method for 
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similar kind of problems and obtained accuracies of order 10~^^°. Problems which consist of 
systems of coupled ODE's with multiple delta function potentials can be solved by, for example, 
the method of shooting to multiple fitting points and implementing the delta functions as 
boundary conditions at those points [17]. Moreover, problems involving moving singularities 
are always are difficult to handle. However a few algorithm have been recently introduced to 
solve these problems (see for example [18]). However, even in this simplest case, the success of 
applying any direct numerical integration method depends on the quality of initial guesses for 
the boundary conditions and energy eigenvalues. Moreover, one usually encounters difficulties 
with the intrinsic instabilities of typical problems, and rarely with the existence of actual 
solutions which posses rapid oscillation. 

Here we discuss an alternative technique for finding the bound states of the time- independent 
Schrodinger equation, applicable for any potential which supports such states. This method 
has the following seven distinct advantages: It is very simple, fast, can be extremely accurate, 
does not have the aforementioned difficulties with the choice of boundary conditions. It is very 
robust and stable, i.e. it does not have the instability problems due to the usual existence 
of divergent solutions of most physical problems. These problems usually produce difficulties 
for the spatial integration routines such as Finite Difference Method (FDM). Moreover, we 
have encountered problems in quantum cosmology whose exact solutions posses very rapid 
oscillations [19] which prevent any successful apphcation of more routines such as FDM, and 
we were able to solve this problem using our alternative technique with ease [20]. Here, our 
main goal is to improve this method by an optimization procedure which first of all allows us 
to define the errors uniquely and secondly decrease them drastically. This method can also 
easily handle cases with mild moving singularities, which also occurred in the aforementioned 
problem. Finally, and perhaps most importantly, we can obtain the wave functions and en- 
ergies of as many of the bound states as desired with a single run of the algorithm. This 
method, which was first introduced by Galerkin over 90 years ago, consists of first choosing 
a complete orthonormal set of eigenstates of a, preferably relevant, hermitian operator to be 
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used as a suitable basis for our solution. For this numerical method we obviously can not 
choose the whole set of the complete basis, as these are usually infinite. Therefore we make 
the approximation of representing the solution by a superposition of only a finite number of 
the basis functions. By substituting this approximate solution into the differential equation, 
a matrix equation is obtained. The energies and expansion coefficients of these approximate 
solutions could be determined by the eigenvalues and eigenfunctions of this matrix, respec- 
tively. This method has been called the Galerkin Method, and is a subset of the more general 
Spectral Method (SM) [21, 22, 23]. Spectral methods generally fall into two broad categories. 
The interpolating, and the noninterpolating method. In the first category, which includes the 
Pseudospectral and the Spectral Element Methods, one divides the configuration space into a 
set of grid points. Then one demands that the differential equation be satisfied exactly at a set 
of points known as the collocation or interpolation points. Presumably, as residual function 
is forced to vanish at an increasingly larger number of discrete points, it will be smaller and 
smaller in the gaps between the collocation points. The noninterpolating category includes the 
Lanczos tau-method and the Galerkin method, mentioned above. The latter is the method 
that we use and, in conformity with the usual nomenclature, we shall simply refer to it as the 
Spectral Method. The interesting characteristic of this method is that it is completely distinct 
from the usual spatial integration routines, such as FDM, which concentrate on spatial points. 
In SM the concentration is on the basis functions and we expect the final numerical solution to 
be approximately independent of the actual basis used. Moreover in this method, the refine- 
ment of the solution is accomplished by choosing a larger set of basis functions, rather than 
choosing more grid points, as in the numerical integration methods. We should note that we 
are implicitly assuming that the true solution is expandable in any complete orthonormal basis 
such as the Fourier basis. However, first of all this requirement is usually satisfied for cases of 
physical applications, secondly we have found that when this requirement is not satisfied the 
method does not fail, but wc loose overall accuracy. 

In this paper we concentrate on the Fourier basis and refine SM by introducing three distinct 
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optimization methods. Then we show that these three procedures, perhaps surprisingly, yield 
identical results which we shall henceforth call the Variationally Improved Spectral Method 
(VISM). This refinement adds two particularly important advantages to SM: first we find a 
unique method to define the error for the cases which are not exactly solvable, and when tested 
on the exactly solvable cases matches closely with the standard definition of the error. Second, 
This refinement, as we shall show, usually improves the accuracy drastically. Then we show 
an application of this amazingly powerful method to the problem of finding the bound states 
of the time independent Schrodinger equation. 

The remainder of this paper is organized as follows. In Section 2, we present the underlying 
theoretical bases for the formulation of the SM in connection with the problems of quantum 
nature. Then we introduce our optimization procedures to obtain VISM. In Section 3, we 
first use this method for the Simple Harmonic Oscillator (SHO) , which is an exactly solvable 
problem for three important reasons: first to illustrate the method, second to show that the 
three distinct optimization procedures yield the same results, and third to test the method 
thoroughly. We then apply this method to two perturbed harmonic oscillators, the first with 
quartic anharmonic term, and second with a rapidly oscillating trigonometric anharmonic 
term. Neither problem is exactly solvable, and are particularly chosen to illustrate some 
of the powerful features of this method. We compare the results for the earlier case with 
those obtained by the usual first order perturbation theory method, the conventional and a 
variationally improved Sturmian approximation method [24], and a highly accurate method 
[25]. The latter, though in principle an approximate method, can determine rather precisely 
the energy levels for this problem from the quantization of an angle variable with an accuracy 
of order 10~^°. The trigonometric anharmonic problem has some very interesting physical and 
mathematical aspects which we shall point out. In Section 4, we state our conclusions and 
some final remarks. 



5 



2 The Spectral Method 

Let us consider the time-independent one-dimensional Schrodinger equation, 

where m, and £^ stand for the reduced mass, potential energy, and energy, respectively. 

Throughout this paper, we only examine the bound states of this problem, i.e. the states 
which are the square integrable. Therefore the general ODE that we want to solve is a linear 
one that can be written in the form. 



where. 



+ f{x)ij{x)=eij{x), (2) 



As mentioned before, any complete orthonormal set can be used for the SM. We use the Fourier 
series basis as an example. That is, since we need to choose a finite subspace of a countably 
infinite basis, we restrict ourselves to the finite region —L < x < L. This means that we can 
expand the solution as, 

H^) = E E 9i (^) , (4) 
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In the above choice of the basis we are implicitly assuming periodic boundary condition. It 

is interesting to note that since we are interested only in the bound states, we could choose 

only the sine series but with the shifted problem x — > x — L. Our shifted domain will be 

< X < 2L and the potential function is shifted accordingly. That is our new basis functions 

and the differential equation are 
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For simplicity we shall refer to this alternate approach the confinement boundary condition. 
Now going back to our original formulation we can also make the following expansion, 

f{x)i.{x)=j:Y.Bm, 9^(^), (7) 
i m \ L J 

where -B^^j are coefficients that can be determined once f{x) is specified. By substituting Eqs. 
(4,7) into Eq. (2) and using the differential equation of the Fourier basis we obtain. 



E 
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j = £ E ^rn,i 9i [-J^j . (8) 



Because of the linear independence of 5'i(^^), every term in the summation must satisfy. 
It only remains to determine the matrix B. Using Eq. (7) and Eq. (4) we have, 

Y.Bm,i9ii-^) =T.^m,im9i(-r)- (10) 

m,i ^ ^ ^ m,i \ L J 

By multiplying both sides of the above equation by 5'i'(^^^) and integrating over the x-space 
and using the orthonormality condition of the basis functions, one finds, 

Bjn,i = E ^rn',i' j 9i f (x) gi> (Ull^j ^^^^ ^ ^ ^i>Cm,m' , (H) 

m',i' ~^ V / m',i' 

Therefore we can rewrite Eq. (9) as, 

~~r' ) ^m,i + 5^ Cm,m',i,i' ^m,' ,i' = ^ ^m,i- (12) 
^ ^ m',i' 

Where the coefficients Cm,m',i,i' are defined by Eq. (11). It is obvious that the presence 
of the operator f{x) in Eq. (2), leads to nonzero coefficients Cm,m',i,i' in Eq. (12), which in 
principle could couple all of the matrix elements of A. Therefore we have to resort to a 
numerical solution. In general the number of basis elements are at least countably infinite. 
The aforementioned coupling of terms in the main matrix Eq. (12) forces us to make the 
approximation of using a finite basis. It is easy to see that the more basis functions we 



include, the closer our solution will be to the exact one. By selecting a finite subset of the 
basis functions, e.g. choosing the first 2N which could be accomplished by letting the index 
m run from 1 to A?" in the summations, equation (12) can be written as, 

DA^sA, (13) 

where D is a square matrix with (2A'") x (2A'") elements. Its elements can be obtained from 
Eq. (12). The eigenvalues and eigenfunctions of the Schrodinger equation are approximately 
equal to the corresponding quantities of the matrix D. That is the solution to this matrix 
equation simultaneously yields 2A'" sought after eigenstates and eigenvalues. The only problem 
which remains is to solve the eigenvalue problem Eq. (13), and to control the round-off errors. 
This is often a serious issue for the usual spatial integration method using double precision. 
However, we can easily overcome this problem and obtain an extremely high precision. This is 
possible if the following conditions are satisfied. First, the potential energy and its derivatives 
should be smooth. Second, the programming language (such as MATHEMATICA, MAPLE, 
etc) should be capable of keeping the desired number of significant digits. Third, when the 
integrations involved in the calculations (Eq. 11) can be done analytically the coefficients in 
the matrix equation (Eq. 12) become exact and the accuracy increases. This also reduces the 
computational time drastically. At this point it is worth mentioning that when the Hamiltonian 
commutes with the parity operator, we can simplify our task of solving Eq. (13) by separating 
the search for the positive and negative parity solutions. 

Now we can introduce our optimization procedure. We are free to adjust two parame- 
ters: 2A'", the number of basis elements used and the length of the spatial region, 2L. This 
length should be preferably larger than spatial spreading of all the sought after wave func- 
tions. However, if 2L is chosen to be too large we loose overall accuracy. As we shall show, 
the error decreases extremely rapidly as the number of basis elements is increased. However, 
it is important to note that for each A^, L has to be properly adjusted. This is in fact an 
optimization problem and is not a trivial task and requires some further analysis. We shall 
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denote this optimal quantity by L{N). We have come up with three methods to determine 
this quantity. These are based on the minimization of any of the eigenvalues (usually chosen 
to be those of the ground state), its error, or the error of the corresponding eigenf unction. For 
the minimization of the energy we proceed as follows. For a few fixed values of we compute 
E[N, L) which invariably has an inflection point in the periodic boundary condition and a 
minimum in the confinement boundary condition. Interestingly enough this points exactly co- 
incide. Therefore, all we have to do is to compute the position of these infiection or minimum 
points and compute an interpolating function for obtaining L{N). Obviously the more points 
we choose the better our results will be. We shall explain the other two optimization methods 
in section 3. Once L{N) is computed, we can use the method described above to compute the 
eigenvalues and eigenfunctions. The necessary program is usually about 15 lines. As we shall 
see, the addition of this refinement to SM can have dramatic consequences. Throughout this 
paper we use VISM. 

Computation of the relative error in the exactly solvable cases is straightforward. For 
example for computing the relative error of the eigenvalue, denoted by Se, we only need to 
find the absolute value of the difference between the result and the exact one and divide by 
the latter. For cases which are not exactly solvable, wc compute the difference between the 
eigenvalues for a given N and those obtained with + 1, both lying on the L{N) curve. We 
shall denote the error computed by this procedure Se- We have computed L(N) for all cases, 
and subsequently computed the eigenfunctions, eigenvalues and their errors using this method, 
and checked their validity in the exactly solvable case of SHO. Obviously to obtain consistent 
results we have to keep the same precision throughout the calculations. 

At this point we should mention that the only weakness of this method that wc have found 
is that, like most other routines, the more discontinuous the potential or its derivatives, the 
less accurate our solution will be. This is due to the fact that these discontinuities would 
induce associated discontinuities in the wave functions or their derivatives via the Schrodinger 
equation. In these cases we would need more basis functions, in particular high frequency ones, 
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to reproduce these features of the wave functions. 

3 Some applications of the Spectral Method 

In this section, for illustrative purposes, we first apply VISM to find the bound states of a 
Simple Harmonic Oscillator (SHO) which is an exactly solvable case, for the reasons stated 
at the end of the Introduction section. Briefly, the purpose is to readily check the vahdity of 
our three distinct optimization procedures, which includes our prescriptions for finding L{N), 
and compare them. We can also check the overall accuracy of our results. We then apply this 
method to two perturbed harmonic oscillators, the first with quartic anharmonic term, and 
second with a rapidly oscillating trigonometric anharmonic term. Neither problem is exactly 
solvable. We compare our results for the earlier case with some other reported results. 

a. Simple Harmonic Oscillator 
The Schrodinger equation for a SHO is, 

^ ' + -rmj^x^il){x) = E^{x), (14) 



2m dx'^ 2 

where u is the natural frequency of the oscillator. Dividing both sides by huj/2, we convert 
this differential equation into the following dimensionless form, 

-f^^ti^ + x'^M)^ E'M), where x' ^ .f^x, E' = ^E. (15) 

This differential equation is exactly solvable and its eigenvalues and eigenf unctions, which are 
all bound states, can be easily found analytically and are well known, 

K = (n + \)hu^, M^) = (-) ^^^^e— n = {0, 1, 2, ...}, (16) 

where Hn{x) denote the Hermite polynomials. Using VISM we can calculate approximately the 
energy levels and the corresponding eigenfunctions of this Hamiltonian. The computation of 
the errors of the wave functions are analogous to that of the energy. We divide the configuration 
space into M grid points. Then, we average the square of the absolute value of the difference 
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between the exact solution and that obtained by the VISM on the grid points, 

_ Inexact _ ^5M| _ E.^il^r^HO-^^"^HOr 

In the left part of figure 1 we show the ground state energy computed using SM for the fixed 
value of the iV = 5 as a function of L using periodic boundary condition. In the right part of 
figure 1 we show the corresponding quantity computed using confinement boundary condition. 
Note the existence of the inflection and minimum points which exactly correspond to each other 
in the figure. These points determine Lib). We repeat this procedure for a few other values 
of N. After plotting these values we can obtain an interpolating function L{N). Our second 
optimization method is based on the minimization of the error of the eigenvalue {Seq), and to 
proceed we first exhibit a semi-log plot of the square of the exact error for the ground state 
energy using SM in terms of N and L in Fig. 2. Note the existence of a valley in this figure 
indicating the optimal quantity L{N), which can gives us the best values for the eigenvalues 
and eigenfunctions using VISM. Our third optimization method is based on the minimization 
of the error of the wave function {S^q), and proceeds analogously to the previous method. In 
Fig. 3 we show our results for L(N) obtained using the three aforementioned methods and 
their interpolating function. Having determined L{N), we can proceed to compute the bound 
states. We have checked the validity of our results for the eigenvalues and their errors, and 
the eigenfunctions using this method as compared to the corresponding exact values. Table 1 
shows the results for the first 10 eigenfunctions for N = 100. Note the outstanding accuracy 
of Se ~ 10~^^° and 5^^ fa 10~^°. The values for 5e are also shown in Table 1. The values for 
could similarly be calculated, but are not shown here. Now we want to exhibit explicitly 
the errors of the ground state wave function, whose value is shown in Table 1 for N = 100, for 
some other values of N. The left part of Figure 4 shows the exact and approximate ground 
state wave functions for N — {3, 5, 7} using SM with fixed un-optimized L — 10. The right 
part of the same figure shows the exact and approximate ground state wave functions for 
N = {1, 2} using VISM with optimized L = {2.52479, 3.04635}, respectively. Note that in the 
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Figure 1: Ground state energy for SHO versus L for = 5, using SM in units where fiuj = 2. 
In the left part the results using periodic boundary condition are shown and on the right using 
confinement boundary condition. 

un-optimized case (SM) for > 7 the exact and approximate wave functions are practically 
indistinguishable on the graph, while in the optimized case (VISM) this occurs for > 2. This 
clearly shows the significance of our variational improvement to SM. Here, we can report that, 
using VISM, Seo = {2 x 10"^, 4 x 10"^ 9 x 10"^ 2 x IQ-^^} for N = {2, 3, 5, 7}, respectively. 
Extrapolating from this, our minuscule errors for N = 100 are easily justified. In the left 
part of Fig. 5 we show a semi-log plot of the error for the ground state energy, obtained 
using VISM, in terms of N. As can be extrapolated from the figure significantly smaller 
errors are easily obtainable. The interesting feature of this graph is that the error decreases 
exactly exponentially and, interestingly enough, the error of the simple Fourier expansion of 
smooth functions has exactly the same behavior. We interpret this as a strong signal that our 
procedure works well. 

In Fig. 5 we show a semi-log plot of the corrected computation time versus N. Note that 
the curvature of this plot is negative and considering this together with the plot of \n{5E) 
shows how the effective efficiency of our program increases with A^. 

b. Anharmonic Oscillator with a quartic term 
Now we apply this method to an anharmonic oscillator which has a quartic term. This is 
probably one of the most famous problem in quantum mechanics which is not exactly solvable 
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Figure 2: Semi- log plot of the square of the error of the ground state energy of SHO (in units 
where Tioj = 2) versus and L using SM. Superimposed on top of the figure is the L{N) 
obtained from this figure, which indicates the direction of the valley. 
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Figure 3: L versus computed by our optimization methods involving minimization of Eq in 
confinement boundary condition (boxes), minimization of 6Eq (circles), as shown in Fig. 2, 
and minimization of 6^^ (diamonds). Also shown in the figure is the interpolating function 
LiN). 
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Table 1: The results for the first 10 eigenvalues and eigenf unctions of the SHO in units where 
huj = 2. E^^^ are values obtained using VISM with = 100. The values of Se and refer 
to the difference between the quantities computed using VISM and the corresponding exact 
quantities as given by Eq. (17). Note the good correspondence between 6e and 6e shows 
the consistency of our method. Obviously we could not display all the digits of E^^ due to 
insignificant errors. The total computation time was about 50 seconds on a typical Pentium 
2.4 GHz machine for obtaining the first 200 eigenvalues and eigenfunctions, which were all 
obtained with a same single run of the algorithm. It is worth mentioning that choosing N=30 
we obtain 6e ~ lO-i^ with computational time of 0.5 Sec. 
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Figure 4: Left, the exact and approximate ground state wave functions of SHO for N = 
{3, 5, 7} using SM with fixed un-optimized L = 10. Note that for > 7 the exact and 
approximate wave functions are almost indistinguishable on the graph. Right, the exact and 
approximate ground state wave functions of SHO for N = {1,2} using VISM with optimized 
L — {2.52479, 3.04635}, respectively. Note that for N > 2 the exact and approximate wave 
functions are practically indistinguishable on the graph. Also note the drastic effect of our 
optimization procedures. 
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Figure 5: Left, semi- log plot of the error for the first eigenvalue of the SHO obtained by VISM 
using various number of basis functions (2A'" equals the number of basis functions used). Right, 
semi-log plot of the computation time versus N for the ground state of SHO. 

and is used as at least as a toy model. The Schrodinger equation for this model is. 



Using VISM we can find the bound states energy spectrum and the corresponding eigenfunc- 
tions of this Hamiltonian. The results that we have obtained using N — 100 are extremely 
accurate {6eo ~ 10"^^°) (Table 2). We have also compared our method and results with some 
other more or less routine methods such as the zero and first order perturbation theory (Table 
2), the conventional Sturmian approximation of Ref. [26], the zero, first and second order 
variational sturmian approximation of Ref [24], and the highly accurate values of Ref. [25]. 
The results of these comparisons is that the seven advantages of our method mentioned in the 
introduction are clearly justified. We just have to mention that the reported accuracy in Ref. 
[25] was only 90 significant digits. The accuracy of the other methods were even lower, that is 
not better than 0(10"^). 

c. Harmonic Oscillator perturbed by rapid oscillations 
Now we want to solve an example which exhibits another less explored and powerful feature 
of VISM. As mentioned in the introduction this method can easily handle problems whose 
solutions exhibit rapid oscillations, in sharp contrast to the spatial integration methods. A 
rather interesting example that we have constructed for this purpose is a harmonic oscillator 
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26.505554752536617417469503006738723676057932189542 


110 


19 


0.28 


32.58 


0.23 



Table 2: The first 10 energy levels of the Schrodinger equation (Eq. (18)) whose dimen- 
sionless form is {—d^jdx''^ + x'"^ + ^ x"^)'4}{x') = E'ip{x'). We have chosen the parameter 
e' — Ae/{muj^) — 1/10 and exhibited the results for the energies in units where hu — 2. 
E^^^ are the values obtained using VISM with = 100, and SD denotes the significant dig- 
its. For space limitations only the first 50 significant digits arc displayed. E^f^ and E^^^^ denote 
the energy eigenvalues obtained using zero and first order perturbation theory, respectively. 
Obviously the accuracy of these methods are far inferior compared to VISM. 



perturbed by rapid oscillations whose precise Schrodinger equation is given by. 



+ {^muj'^x'^ -\- acos{(5'Kx)^ ijj{x) — Eijj{x), 



(19) 



2m dx"^ 

where cu is the natural frequency of the oscillator and a and /3 are arbitrary constants. This 
differential equation is not exactly solvable and for high frequency perturbations j3 the spatial 
integration routines like FDM would have serious difficulties in these situations. They should 
consider many spatial points for overcoming the rapid oscillations of the potential which in- 
creases the time and round off errors of the routine and decreases the efficiency and stabihty 
of the method. However, this is not a case for the VISM which can also handle these type of 
potentials very easily. For large /3 the behavior of the potential is very oscillatory and centered 
around the curve ^miux^. As expected all the eigenstates of this problem are bound states. 
The results for the ground state are shown in the Fig. 6. In the left part of the figure we show 
the full potential U{x), the ground state wave function {il)o{x)), and a zoomed box highlighting 
the fine structural behavior of the wave function. In the right part of the figure we show the 
ground state energy Eq versus A^. Note for TV smaller than (100 here) VISM is not sensi- 
tive enough to respond to the rapidly oscillating part of the potential and the results are very 
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Figure 6: Left, the potential of the harmonic oscillator perturbed by rapid oscillations, whose 
Schrodinger equation (Eq. (19)) in dimensionless form is {—d'^/dx''^+x'^+a'cos{f3'nx'))ilj{x') = 

E'ip{x'). We have chosen the parameters P' = ^j2/m(3/uj = 10, and a' = {2/hu)a = 10 and 
exhibited the results for the energies in units where hcu — 2. Superimposed on the same graph 
is the ground state wave function calculated using VISM with = 150. The zoomed box 
exhibits the behavior of the wave function and the potential for —0.7 < x < 0.7. Right, the 
ground state energy versus N for the same parameters as in the left figure. 

close to those of the (unperturbed) SHO. As is apparent from the figure, for N even slightly 
larger than /?L, the VISM easily incorporates the rapidly oscillating part of the potential and 
rapidly approaches the exact energy eigenvalue as N increases. This explains the sharp drop 
in the graph of Eq(N). In particular for N = 150 we obtain ^g,, ^ 10^^°. The physics of 
this phenomenon is very clear: When the set of basis functions is large enough to include 
those whose frequencies are at least as large as the ones induced in the wave function by the 
oscillations or discontinuities in the potential, via the Schrodinger equation, VISM can easily 
finds the accurate eigenfunctions and the corresponding eigenvalues, the fine structure of the 
ground state energy exhibited in the left part of the figure is also very interesting. The shght 
ripples in the wave function are produced by the conflicting tendencies of the wave function 
to accumulate in the rapidly occurring wells of the perturbed part of the potential and that of 
the second derivative operator to smooth out the wave function. 
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4 Conclusions 



We have introduced three distinct yet equivalent variational improvements for the Fourier 
Spectral Method and used them as an extremely accurate method for obtaining the energies 
and wave functions of the bound states of the time-independent Schrodinger equation. In this 
method a finite basis is used for approximating the solutions. Our variational improvement of 
the method is to calculate an optimized spatial domain for a given number of basis elements, 
denoted by L,{N). The three optimization procedures are based on minimizing the eigenvalues, 
the error of the eigenvalues, and on minimizing the error of the eigenstates. Our refinement 
schemes usually improve the accuracy of SM drastically, as can be seen in Fig. 4, for example. 
This improvement increases rapidly with N. In particular, when the problem is not exactly 
solvable there does not seem to be any other canonical way to fix L and the errors. We 
applied this method to an exactly solvable problem and easily found an extraordinarily good 
agreement with the exact solutions (errors of order 10~^^°). In the quartic anharmonic oscillator 
case which is not exactly solvable, the accuracy of this method is much higher than some of 
the more conventional methods such as the perturbation method, the conventional sturmian 
approximation, and the variational sturmian approximation. In the problem of SHO perturbed 
by trigonometric anharmonic term, we observed how easily and accurately VISM handles 
problems involving potentials with very rapid oscillations. To summarize, this method is 
very simple, fast, extremely accurate in most cases, very robust and stable, can easily handle 
solutions with rapid oscillations and moving singularities, and there is no need to specify the 
boundary conditions on the slopes. Most importantly, one can obtain the energies and the 
wave functions of as many of the bound states as desired with a single run of the algorithm. 
The main sources of error are using too few number of basis elements and having potentials 
with major discontinuities. When the latter dose not exist we can obtain extraordinary good 
results, e.g. 130 significant digits, by choosing appropriate parameters. This method can be 
easily extended to D dimensional Schrodinger equation. For example. We have solved 2-D 
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Harmonic Oscillator V{x,y) = x'^ + y"^ and QCD potential V{x,y) ~ x^y^ using VISM and 
found extremely highly accurate results [27]. In 2-D SHO case the accuracy with 15 significant 
digits can be easily obtained with only 22 basis functions, and this shows the efficiency of 
VISM. 
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